"""
This file contains the code used for obtaining Theorem 17.

Main method: 
    short_fillings_with_long_systole(m_name, m = None, length = 7.515, min_systole = 0.163 , verified = True, bits_prec = 202, verbose = False)

Input: m_name, the string of the name of a 1-cusped manifold in the census (e.g. m_name = 'o10_000001')
Output: a list of pairs (p,q) representing the slopes containing all slopes that are
    (1) no longer than length (or indistinguishable under the given precision)
    (2) m(p,q) has systole length at least min_systole (this is rigorous)

Detailed logs will be printed when setting verbose to True
"""

exceptional_path = '../exceptional_census/exceptional_fillings.dat'
# the path to the file saving the exceptional fillings,
# whose relative path in this repository is set by default as above.

trys = 200
# for some examples, it is necessary to set trys = 3000 to make the computation successful.

import ast
from snappy_10_tets import snappy
# version of snappy used: 3.2
# version of snappy_10_tets used: 1.20

er_count = 0
er_limit = 30

excep_slope_dict = dict()

with open(exceptional_path, 'r') as f:
    for line in f:
        line = line.strip()
        [m_name, slope] = line.split('(')
        slope = ast.literal_eval('('+ slope)

        if m_name not in excep_slope_dict.keys():
            excep_slope_dict.update({m_name:[slope]})
        else:
            excep_slope_dict[m_name].append(slope)

hard_set = set(['m007(3,1)', 'm135(1,3)'])
hard_slope = set(['o10_107134(1, 0)'])
# Some examples encountered whose geometric triangulations are hard to find,
# making verified computations hard and require other methods

def closed_isosigs(snappy_manifold, trys=trys, max_tets=60, verbose = False):
    """
    Generate a slew of 1-vertex triangulations of a closed manifold
    using SnapPy.
    
    >>> M = snappy.Manifold('m004(1,2)')
    >>> len(closed_isosigs(M, trys=5)) > 0
    True
    """

    M = snappy.Manifold(snappy_manifold)
    assert M.cusp_info('complete?') == [False]
    surgery_descriptions = [M.copy()]

    try:
        for curve in M.dual_curves():
            N = M.drill(curve)
            N.dehn_fill((1,0), 1)
            surgery_descriptions.append(N.filled_triangulation([0]))
    except snappy.SnapPeaFatalError:
        pass

    if len(surgery_descriptions) == 1:
        # Try again, but unfill the cusp first to try to find more
        # dual curves.
        try:
            filling = M.cusp_info(0).filling
            N = M.copy()
            N.dehn_fill((0, 0), 0)
            N.randomize()
            for curve in N.dual_curves():
                D = N.drill(curve)
                D.dehn_fill([filling, (1,0)])
                surgery_descriptions.append(D.filled_triangulation([0]))
        except snappy.SnapPeaFatalError:
            pass

    ans = set()
    for N in surgery_descriptions:
        for i in range(trys):
            T = N.filled_triangulation()
            if T._num_fake_cusps() == 1:
                n = T.num_tetrahedra()
                if n <= max_tets:
                    ans.add((n, T.triangulation_isosig(decorated=False)))
            N.randomize()
            if verbose:
                print('\033[K' + f'closed_isosigs, Trial: {i+1}, Ans number: {len(ans)}', end = '\r')

    if verbose:
        print()
    return [iso for n, iso in sorted(ans)]

def systole_length(m, verified = True, bits_prec = 202, verbose = False):
    global er_count
    global er_limit

    er = ''

    try:
        sigs = closed_isosigs(m, verbose= verbose)
        assert len(sigs) >0

        if verbose:
            print(f'closed_isosigs obtained: {len(sigs)} triangulations')

        i = 0 
        for sig in sigs:
            i += 1
            if verbose:
                print('\033[K' + f'Try computing systole length with the {i}th triangulation', end = '\r')
            n = snappy.Manifold(sig)
            n.simplify()
            try:
                ans = n.length_spectrum_alt(count = 1, verified = verified, bits_prec = bits_prec)[0].length.real()
                num_tet = snappy.Triangulation(sig, remove_finite_vertices = False).num_tetrahedra()
                if verbose:
                    print(f'Systole length successfully computed with the {i}th triangulation {sig} with {num_tet} tetrahedra')
                return ans
            except Exception as e:
                er = e
                pass

        raise er
    except Exception as e:
        if verbose:
            print()
            print(f'Computation failed, retrying systole_length. er_count = {er_count}')
            print(f'Error message: {e}')
        er_count += 1
        if er_count< er_limit:
            return systole_length(m,verified = verified, bits_prec =  bits_prec, verbose = verbose) 
        else:
            raise RuntimeError 


def short_fillings_with_long_systole(m_name, m = None, length = 7.515, min_systole = 0.163 , verified = True, bits_prec = 202, verbose = False):
    global excep_slope_dict
    global hard_set
    global er_count
    
    if m is None:
        m = snappy.Manifold(m_name)

    length = RealIntervalField(bits_prec)(QQ(length))
    min_systole = RealIntervalField(bits_prec)(QQ(min_systole))

    try:
        area = m.cusp_areas(verified = verified, bits_prec = bits_prec)[0]

        slope_list = m.short_slopes(length = length*area, verified = verified, bits_prec = bits_prec)[0]
        # this might give slopes longer than 7.515 (normalized), but we filter the list rigorously below by systole length

        ans_list = []
        if verbose:
            print(f'A total of {len(slope_list)} short slopes obtained for {m_name}')
        
        for slope in slope_list:
            er_count = 0
            if verbose:
                print(f'Examinating systole of {m_name}{slope}. Skip if filling is exceptional or hard.')
            if (m_name not in excep_slope_dict.keys() or slope not in excep_slope_dict[m_name]):
                if (m_name + str(slope)) in hard_slope:
                    if verbose:
                        print(f'Adding {slope} into answer list')
                    ans_list += [slope]
                    continue
                identified = snappy.Manifold(m_name + str(slope)).identify()
                if identified and str(identified[0]) in hard_set:
                #if m_name in hard_slope_dict.keys() and slope in hard_slope_dict[m_name]:
                    if verbose:
                        print(f'{m_name}{slope} identified to be one of the hard manifolds, adding slope to list')
                    ans_list += [slope]
                else:
                    #print('check slope ' + str(slope))
                    n = snappy.Manifold(m_name + str(slope))
                    if n.solution_type() == 'all tetrahedra positively oriented': 
                        # compute only with geometric triangulation to avoid mistakes
                        try:
                            if verbose:
                                print('Checking core_length')
                            core_length = n.cusp_info('core_length')[0].real()

                            if len(str(core_length)) > 6 and RealIntervalField(bits_prec)(QQ(core_length)) < min_systole:
                                if verbose:
                                    print(f'{m_name}{slope} has core length {core_length}, shorter than desired. Proceed to the next slope.')
                                continue

                            if verbose:
                                print('Checking dual curves')
                            n_curves = n.dual_curves()

                            if len(n_curves) > 0 and len(str(n_curves[0]['filled_length'].real())) > 6 and RealIntervalField(bits_prec)(QQ(n_curves[0]['filled_length'].real())) < min_systole:
                                if verbose:
                                    print(f'{m_name}{slope} has dual curve length {n_curves[0]['filled_length'].real()}, shorter than desired. Proceed to the next slope.')
                                continue
                        except:
                            pass
                    
                    if verbose:
                        print(f'Computing the systole length of {m_name}{slope}')
                    this_systole = systole_length(n, verified = verified, bits_prec = bits_prec, verbose=verbose)

                    if not this_systole < min_systole:
                        if not this_systole > min_systole:
                            if verbose:
                                print(f'Cannot decide if the systole of {m_name + str(slope)} is longer than desired, asserting False and retry from draft')
                            assert False
                        # using assert, if the above is not the case,
                        # we discard the current answer and redo the computation from draft.
                        
                        if verbose:
                            print(f'Adding {slope} into answer list')
                        ans_list += [slope]

        return ans_list
    except Exception as e:
        if verbose:
            print(f'Computation failed, retrying from draft. er_count = {er_count}')
            print(f'Error message: {e}')
        er_count += 1
        if er_count< er_limit:
            m.randomize()
            return short_fillings_with_long_systole(m_name, m, length, min_systole, verified, bits_prec, verbose=verbose)
        else:
            raise RuntimeError